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We study the influence of multiple bands on the properties of Josephson junctions. In particular 
we focus on the two gap superconductor magnesium diboride. We present a formalism to describe 
tunneling at a point contact between two MgB2 electrodes generalizing the transfer-matrix approach 
to multiple bands. A simple model is presented to determine the effective hopping amplitudes 
between the different energy bands as a function of the misorientation angle of the electrodes. We 
calculate the critical current and the current- volt age characteristics for N-I-S and S-I-S contacts with 
different orientation for junctions with both high and low transparency. We find that interband 
tunneling processes become increasingly important with increasing misorientation angle. This is 
reflected in certain features in the differential tunneling conductance in both the tunneling limit as 
well as for multiple Andreev reflections. 
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I. INTRODUCTION 

The two gap superconductor MgB 2 is considered to 
be a good candidate for superconducting devices due to 
its comparatively high critical temperature, its easy han- 
dling and cheap preparation. Presently, its use in Joseph- 
son junctions is being investigatedii2i2i2i£. One of its pe- 
culiar properties is the presence of two distinct super- 
conducting gaps existing on different bands at the Fermi 
surface. These bands are n and a bands arising from the 
Boron p z and p x /p v orbitals, respectively^. These two 
types of bands possess different parity with respect to 
reflection at the Boron plane. This different parity sup- 
presses scattering and transitions between the two types 
of bands, which is thought to be the main reason for 
the exceptional stability of the two gaps against impu- 
rity scattering^. 

In the vicinity of a Josephson junction this parity can 
be broken, if the two MgB2 electrodes are grown with 
a misorientation angle of the crystal c-axis directions on 
both sides of the junction. The size of the misorientation 
angle would provide a means of tuning the strength of in- 
terband transitions between the two sides of the junction. 
This is what we wish to explore in the present work. In 
particular, we want to study the differential conductance 
of such a junction and demonstrate how these interband 
processes show up in the tunneling conductance and in 
Andreev reflections. 

Due to the importance of Josephson junctions in nu- 
merous applications the Josephson effect has been sub- 
ject of intense theoretical study. The current-phase rela- 
tion as well as the oscillation of the electron transport in 
a voltage biased point contact known as the ac- Josephson 
effect has been theoretically investigated for different sit- 
uationaSiiSiii including unconventional superconductiv- 
itjiiS, spin active magnetic barriersi2ii!ii£ and the effect 
of pair breaking due to disorder. It has been shown 
that within the quasiclassical theory interfaces can be 
described by effective boundary conditions^ that can 



be included by mainly two different approaches: the 
scattering-matrix approach and the transfer-matrix ap- 
proachl^. Both are formally equivalent but the range of 
applicability is different and the problem under consid- 
eration can help to decide which one should be chosen. 
Here, we are going to generalize the transfer-matrix ap- 
proach to the case of a multiband Josephson junction and 
use it to calculate the differential conductance of a MgB2 
Josephson junction with misoriented electrodes. The ef- 
fects of quasiparticlc and Josephson tunneling between 
two multiband electrodes with application to MgB2 have 
been theoretically discussed by Brinkman et al. within an 
extended BTK-model 18 explaining convincingly the ab- 
sence of a double gap structure in c-axis tunneling spec- 
tra. Novel Josephson effects between multiband and sin- 
gleband superconductors including the claim of a phase 
shift of 7r between the different gaps on the multiband 
side have been theoretically discussed by Agterberg et 
alJ2 showing the rich and interesting physical content 
of Josephson effects in multiband systems. The exper- 
imental observation of subharmonic gap structures in 
MgB 2 single crystal junctions^ and MgB 2 /Nb micro- 
constrictionsSi due to multiple Andreev reflections can 
be understood within a multiband transfer-matrix ap- 
proach. We will start in the next section with a descrip- 
tion of the transfer-matrix model generalizing the ap- 
proach to include the effect of multiple gaps. In the third 
section we will derive the current expression for both the 
equilibrium case without an applied voltage and the non- 
equilibrium case with an applied voltage. In the fourth 
section we will model the effective hopping amplitudes 
at the interface making use of some general considera- 
tions. In the fifth section we will show the results of our 
calculations while in the last section we will conclude. 



2 



II. THE TRANSFER-MATRIX APPROACH 
FOR MULTIPLE BANDS 



The quasiclassical theory of superconductivity has 
been proven to be a powerful tool to consider spa- 
tially inhomogeneous problems in equilibrium and non- 
equilibrium, for example to calculate the local quasipar- 
ticle density of states around vortices or in the vicinity 
of boundaries. However it is only valid within the qua- 
siclassical limit kp^ ^> 1 and can therefore only describe 
situations with slowly varying fields and potentials (on 
the scale of the coherence length). The interface be- 
tween two superconducting regions of different orienta- 
tion or between a normal metal and a superconductor 
represents a strong pertubation on a lcngthscalc much 
smaller than the coherence length which is - in principle 
- out of the range of validity of the quasiclassical theory. 
To describe this strong pertubation within the quasiclas- 
sical limit one has to find effective boundary conditions 
that connect the solutions on both sides. The first for- 
mulation of these boundary conditions for nonmagnetic 
interfaces has been found by Zaitse\iit Afterwards they 
have been generalized by Millis et ali2£ for magnetic ac- 
tive interfaces and have been explicitly solved by She- 
lankov— for equilibrium problems and complemented by 
EschrigSi for non-equilibrium problems within the power- 
ful Riccati parametrization of the quasiclassical theory 25 . 
To describe the complicated processes at the interface be- 
tween two superconductors including multiple Andreev 
scattering the transfer-matrix approach has been proven 
to be as well suitable as the scattering-matrix approach 
but gives a more intuitive understanding. To describe 
what is happening at the boundary between two super- 
conductors with several Fermi surfaces we will generalize 
in this work the transfer-matrix description of Cuevas 
and Fogelstrom 12 and Kopu et almS for multiple Fermi 
surfaces. The basic idea of this approach is to consider 
two decoupled superconducting electrodes described by 
two electrode Green's functions goo,L/R. that are calcu- 
lated in a left (L) and a right (R) half-space with an 
inpenetrable surface. In a second step we allow virtual 
hopping processes between the two electrodes including 
the possibility of intraband and interband hopping. This 
serial expansion is exactly summed up to infinite order 
using the technique known from the Dyson formalism. 
As has been shown by Kopu et al. this pertubative 
transfer-matrix approach is equivalent to other methods 
and can be used to calculate interfaces with arbitrary 
transparency. In the following we will repeat the steps 
described by Cuevas and Fogelstrom and we will discuss 
the consequences that arise from the multiband character 
of the superconducting pairing interaction. To describe 
the hopping processes we start with a phenomenological 
transfer Hamiltonian 



= E E iftJ + ( £ r,S<i^C (i) 



where a and a' denote band indices corresponding to a 
special Fermi surface and introducing the effective hop- 
ping amplitudes w£^, describing hopping processes from 
band a' in the right electrode to band a in the left elec- 
trode, and UjI" correspondingly. In this formulation we 
have enlarged our Hilbert space to include not only the 
band index but also the quantum number of the elec- 
trodes. The full Green's function and the T-matrix can 
be written as 
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taking into account reflection (LL and RR) and trans- 
mission processes (LR and RL). The Green's functions 
Gij and T-matrices Tij themselves are N x iV-matrices in 
the enlarged band space where N is the number of bands 
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consisting of 2x2-matrices in Keldysh space, denoted by 
the check symbol 



ij 



6 G% aa 



(4) 



The retarded, advanced and Keldysh Green's functions 
are 2x2 Nambu-Gor'kov matrices in particle-hole space 
while we will neglect the spin degrees of freedom since we 
are concerned with spin singlet systems and we assume 
only non- magnetic interactions. Besides the full Green's 
functions we also have to consider the Green's functions 
of the unperturbed left and right electrodes, which take a 
diagonal form in the space of the two electrode quantum 
numbers, and we can write down the coupling parameters 
as off-diagonal elements in the same notation: 



Goo = 



Gco,L 
Goo,R 



V L R 

vrl 



(•5) 



Since we want to neglect direct interband scattering in 
the decoupled left and right electrodes the self-energy 
and therefore also the unpertubed Green's functions take 
diagonal form in band space 
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(6) 

We will see in the following that the full matrix struc- 
ture of the Vij in band space prevents us from an easy 
decoupling of the different matrix elements as it can be 
done in the space of the electrode quantum numbers. 
Therefore we have to solve the full matrix problem in 
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band space. With regard to the application to the two- 
band superconductor MgB 2 , where interband scattering 
processes are suppressed^, we will go on with these sim- 
plified diagonal electrode Green's functions but of course 
it is straight forward to extend the considerations to a 
full matrix problem including interband scattering. The 
transmission processes can now be summed up in form 
of a T-matrix equation with 



7" = v + v o oT = v + voQov 



(7) 



This equation has to be complemented by the Dyson 
equation which reads 

Q = Goo + Goo ° T ° Goo i = Qoo + Goo ° v o g (8) 

Using the first identity of Eq. (JJJ repeatedly we can write 



T = V + V O t/oo O V + V O C/qo O V O C/oo o T 



(9) 



If we now employ the diagonal structure of the unper- 
turbed Green's functions und the off-diagonal structure 
of the hopping matrices in the space of the electrode 
quantum numbers we get the following set of decoupled 

equations for the components of T. We note the diagonal 
(reflection) elements of the T-matrix as 

TLL = V L R O Goo.R ° vrl 

+VLR O Goo.R O Vr L o Q^x o f LL (10) 
Trr = Vrl O Qoo.l o V L R 

+vrl Goo,l°vlr°Goo,r T'rr (11) 

and we can write the off-diagonal (transmission) elements 
as 

Tlr = V L R + V L R ° Goo,R ° Vrl ° Goo.L ° T L R (12) 
TRL = Vrl + V RL ° Goo,L ° V L R ° Goo,R ° f RL (13) 

Now we can perform the quasiclassical integration since 
there are no functions with spatial arguments from both 
sides of the interface any more. The quasiclassical £- 
integration reads 

Soo,i(pF,t,t') = ~ J d^fas/pF Goo,i(PF,t,t')y/pF 

(14) 

where the matrix \fpF gives the square root of the density 
of states at the Fermi level in the normal state. This 
matrix is proportional to the identity matrix in Keldysh 
space but has different elements on the diagonal in band 
space 
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We can go on defining the quasiclassical transfer and hop- 
ping matrices as 



i(PF,p' F ) = KyffeVij(j)F,tf F )y/pFfz, (16) 



tij(j>F,PF> t > t ') = ^VPF^ij(PF^F^^)ypF f 30-7) 

Finally we can write down a set of equations defining the 
quasiclassical transfer matrices 
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where the Fermi surface average (■■ - )p F has to be per- 
formed with respect to each different Fermi surface sheet 
under consideration. If we now take a look at the qua- 
siclassical transport equation for each band a where the 
strong perturbation at the interface enters as a source 
term on the right hand side 

i^StipF) + HPf) ~ &i{pF),9i(pF)]l 

= 2nv F,n [iii(^F,PF),goo,i(pF)]'^S[f-f e ) (22) 

complemented by the normalization condition 

g(g>g = -ir 2 i (23) 

we can calculate the rapid change of the quasiclassical 
propagator at the interface by integrating over the source 
term in Eq. (|22H neglecting all the slowly varying fields 
on the left hand side. We get for the difference of the 
quasiclassical Green's function on the incoming and the 
outgoing trajectory - evaluated directly at the boundary 



[&,+ (Pf) -9i,-(PF)] a 

= -2wi [iaipF ,PF),9oo,i(pF)]l 



(24) 



Now it is possible to construct a solution for g^ and 
cji :+ that fulfills the boundary condition given by Eq. I|24l) 
and that also satisfies the normalization condition of the 
quasiclassical propagator given by Eq. 

Si- = 9oo,i + {Sao,i + wrl) ® ki ® (<7oo,; - wri) (25) 
9i.+ = 9oo,i + {goo,i - wri) <8> ki <8> (goo.i + wrl) (26) 

where we have omitted the angular dependencies. Tak- 
ing Si,± as boundary values we can now construct the 
full quasiclassical propagator on incoming and outgoing 
trajectories in the vicinity of an interface with arbitrary 
transparency. In the following we will use it to calculate 



4 



the quasiparticle density of states at the boundary, the 
Josephson currents across the boundary and the currents 
at an applied voltage for a multiband s-wave supercon- 
ductor. Since there are no Andreev bound states at the 
interface without a sign change of the order parameter for 
reflected trajectories we do not expect dramatic changes 
of the pairing potential amplitude near the boundary. In 
this case we can go on without a self-consistent calcula- 
tion of the local pairing potential. 



To proceed further we must take into account that the 
different bands posses different Fermi surface geometries. 
As we have shown in a previous work the three dimen- 
sional tubular network of the tt bands can be approxi- 
matly described by a half torus with a fraction of the 
two toroidal radii of v = 4 while the nearly two dimen- 
sional isotropic Fermi surfaces of the a bands can be de- 
scribed by distorted cylinderg 2 ^^. Taking the interface 
to be orientated perpendicular to the a6-plane direction 
we can write the Fermi surface averages as 



III. CURRENT WITH AND WITHOUT AN 
APPLIED VOLTAGE 

The current density can be calculated integrating the 
Keldysh component of the quasiclassical propagator over 
the quasiparticle energy variable e. We can write the 
current component in the band a as 



3 a (r,t) = eN { F a) 



de 
8ni 



Tr(v F f 3 g K ' a (f,e,t)) a (27) 



In the following we want to calculate only the current 
component across the surface on one side of the contact, 
e.g. the left side, so we can write 
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If we determine the values of - 



4.4- 10 5 s and 



8.2-10 5 — from band structure calculations 18 27 and insert 

s 

the value of v = 4 we finally get the weighting factors for 
the two different bands as 



Dp COS < 



= 1.4- 10° 



,m 



(33) 



where (■ ■ ■} a ,+ means an average on the Fermi surface 
corresponding to the band a with the restriction that 
only angles corresponding to outgoing trajectories with 
V Fn > are taken into account. As we have seen before 
we can calculate the difference of the incoming and the 
outgoing propagator in the vicinity of the surface as 



and 



Tr |f 3 [i LL ,g 00tL \^ ,a \ 



(29) 



This expression can be further simplified in the equilib- 
rium case where the Keldysh propagator can be written 
as difference of the retarded and advanced Green's func- 
tion multipled by a temperature dependent distribution 
function. In the non-equilibrium case where we apply 
a current over the contact we have to calculate the time 
evolution of the transfer-matrix and get a time-dependent 
solution for the current. 



A. Josephson currents 



(^cos</>cos6») T+ = 1.588- 10 5 



(34) 



Now we can calculate the Josephson current as a func- 
tion of the phase difference between the gaps on the left 
and on the right hand side of the junction from Eq. 1281) 
taking the incoming and outgoing Green's functions from 
Eq. 1251) and Eq. (|26|l and evaluating the quasiclassical 
t- matrix as defined in Eq. (|18|l . Since the Josephson cur- 
rent is an equilibrium property the £g>-products reduce 
to simple matrix multiplications. To compare the qua- 
siclassical results of the temperature dependence of the 
critical currents with a simplified model we can write 
down a two band analogon to the Ambegaokar-Baratoff 
formula. Here we have to consider the weighted sum of 
two single band Ambegaokar-Baratoff formulas using the 
temperature dependence of the pairing potentials calcu- 
lated selfconsistently from the two band gap equation: 



To calculate the currents we have to perform the Fermi 
surface average as defined in the previous section. Since 
within our model none of the Green's functions posses an 
angular dependence we only have to evaluate the follow- 
ing expression for the a and for the ir band 



(30) 



Jc (T) = £ eNf ] (v F J a + V«nA^ tanh 



A(«) 
2T 



(35) 

Here T> a is an effective band dependent transmission co- 
efficient. For low transparencies it can be calculated from 

the hopping amplitudes as V a — Air 2 v a ' a 
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B. S-I-S- and iV-Z-S-contacts 

If we apply a constant voltage across an S-I-S Joseph- 
son junction we will find that the phase difference in- 
creases linear according to the Josephson relation <fi(t) = 
<f>Q + ujot where the Josephson frequency is proportional 
to the applied voltage ujq = 2eV/h. This results in an 
oscillating current across the junction and is known as 
the ac- Josephson effect. In this case we have to solve 
the non-equilibrium time-dependent t-matrix equation. 
First of all we will simplify the Keldysh component of 
[III , <?oo,lL and we will show that the problem can 
be reduced to calculate only the retarded and advanced 
parts of the transmission components of the quasiclas- 
sical i-matrix i^ R A and t R j^ '. Then we will determine 
the time dependence of the transfer-matrix components 
mentioned above by expanding them into a Fourier se- 
ries. As shown in the appendix the commutator can be 
transformed in such a way that we can write the spectral 
current density as 
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(36) 



Now we have to find the time evolution of the transfer- 
matrix. Since there is no further complication from the 
multiband character of our problem we can directly follow 
the steps described by Cuevas and FogelstromiS. Start- 
ing with the transfer-matrix equation 



*£ij A (*>*') = v LR + j dti J dt 2 v LR g^x(t,ti) 

^RLfJ^tiM^n^t') (37) 
and Fourier transforming it to energy space 



~R/A, 



O = ^ / & / de'e-^'t^ie, e>) ( :■}* i 



leads to an algebraic equation to determine i_L_R(e, e'). 
Further one can show that the t-matrix allows a Fourier 
expansion of the form 

££V e') = ^fe e + n eV)5{e - e' + neV) (39) 



So the problem of calculating the current can be reduced 
to the problem of calculating the Fourier components 
i nm for the retarded and advanced part of the transfer- 
matrix, that are equally defined as 



inm = tLifi 6 + neV, e + meV) 



(40) 



In the following we will use the same notation for the 
advanced and the retarded part of the i-matrix since they 
are described by the same equation. We will also use 



the notation v = vl R and = v R l. But we have to 
keep in mind, that all t nm are not only 2x2-matrices 
in Nambu-Gor'kov space but are also N x iV-matrices 
in band space. Finally we can write down the algebraic 
i-matrix equation asi£ 

+Vn,n-2tn-2 : m + Vn : n+2t n +2 : m (41) 

where we have defined the matrix components as 



and also 

V n ,n+2 = -vfR,n+\V 

and finally 



fL,n+2 9L,n+2 



Vn.n-2 



ft t ( 

V 9L,n- 







2 JL,n-2 



(42) 



(43) 



(44) 



Here and in the following the notation g in stands for the 
diagonal and the functions fc n and stand for the 
off-diagonal parts of the equilibrium Green's function, 
evaluated at e + neV as 



9i,n = 9ooA e + neV ) 



(45) 



The matrices v nm are only evaluated for n = m ± I and 
we can write them as 



I'm- l,i 
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(46) 



Now the algebraic set of equations defined by Eq. i|4I|) 
can be solved by standard recursive techniques that we 
will briefly demonstrate. A significant simplification can 
be achieved if we use the following identity 



*n,m(e) = *n-m,o(e + meV) 



(47) 



that follows directly from the definition of the t n , m (e). 
Now we can proceed further by defininig some ladder 
operators z^ as 



n,0 



(48) 



It can be easily verified that to,o = 6 is a solution of 
Eq. (|4I|I for n — and therefore all i n fl with even n van- 
ish identically. So we only have to calculate the odd lad- 
der operators. Inserting the above relations in Eq. (|41|) 
leads to recursion formulas for the z^ for n > I or rt < — f 
respectively: 



1 + V„^„±2^„±2 ) V n! „ T 2 (49) 



Since the are getting constant for high values of n 
we can start with z+ =0 and zZ n = counting 
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down to z 3 or up to z_ 3 . After solving the equations for 
t±\fi we are now able to determine all t Kj o by recursive 
application of the ladder operators. To find the t±i,o we 
have to solve two coupled 6x6 matrix equations in the 
combined band and Nambu-Gorkov space: 

h,o = "1,0 + £1,1*1,0 + Vi,-i*-i,o + Vi t3 z^ i-Lfl (50) 

t—1,0 = ^— 1,0 + £—1, -1*-1,0 + V-l,-32_3f-l,0 + V-l,lil,0 

Of course this procedure has to be applied to the retarded 
as well as to the advanced components of the transfer 
matrices. The missing transfer matrices can be found 
using the general relation i^ L R nm (e) = fsi^fmnC 6 )^ as 
for example described by Cuevas and Fogelstrom 12 . With 
this procedure the spectral weight of the non-equilibrium 
current across an S-I-S contact for a given voltage V and 
a given quasiparticle energy e can be calculated. From 
this we get the time dependent non-equilibrium current 
in a voltage biased contact as 

oo 

J« L (V,t)=eN ( F a) (v F ') a , + £ j m e 2mieVt ^ (51) 

m— — oo 

with the Fourier components 

— ^ Tr {f 3 ($LRg%oi^ t0n gL,jLR,nmVLR 
n 

~-lfR ~K fA -A - /r \ 

~ V RL Z RL,0n9L,n l LR,nm9R,m v RL \04) 

, fR ~K TA ? A _ h R fR ~K fA \ 1 

I L LR,0n9R,n l RL,nm9L,m 9 L,0 L LR,0n9 R,n L RL,nm ) ) 

For the special case of an N-I-S contact we can choose 
without loss of generality the right hand side as the nor- 
mal conducting side with = and therefore also 
V n ,n±2 = 0. In this case the ladder operators vanish 
and we only have to calculate the transfer matrix com- 
ponents i±i.o- This can be understood since the lad- 
der operators create repeated Andreev reflections be- 
tween the superconducting electrodes and in the case of a 
superconductor-normal conductor interface we only have 
single Andreev processes. 

IV. DETERMINATION OF THE EFFECTIVE 
HOPPING AMPLITUDES FOR MGB 2 

One of the crucial features of the two gap supercon- 
ductor MgB2 is the orthogonality of the wavefunctions 
belonging to the 7r and a bands, respectively. Here, we 
want to consider the case that the two MgB2 electrodes 
possess a certain misorientation angle of the crystal c-axis 
directions on both sides of the junction. The misorien- 
tation angle allows to change the interband transmission 
between the two sides. In the following we will introduce 
a simple model to determine the effective hopping ampli- 
tudes at the interface between these two misaligned re- 
gions. The a bands are formed by the in-plane a bonds of 




FIG. 1: (color online). The effective hopping amplitudes be- 
tween the different bands of MgB 2 can be modelled by cal- 
culating the overlap integrals of the p x ± ip^-orbitals and the 
Pz-orbitals of the boron atoms, which are responsible for the 
formation of the superconducting a and n bands of this ma- 
terial. The angle 7 measures the relative orientation of the 
two crystals on the different sides of the interface. 



the p x ±ip y orbitals of the boron atoms, while the ir band 
is formed by the weaker tt bonds of the boron p z orbitals 
connecting the different boron layers. To model an in- 
terface between two electrodes with different orientation 
we have to calculate the overlap integrals between the 
different orbitals rotated against each other (see Fig. 
Since all orbitals show rotational symmetry with respect 
to the z-axis it is sufficient to study only the rotation 
around one of the other axes, e.g. the x-axis. We can 
write down the effective one-particle wave function for 
the p x , the p y and the p z orbitals, neglecting the radial 
part of the wave function, as follows 

\Px) = --^(Yi,i(M)-Yi,-i(M)) (53) 

\Py) - -^(YliCMHIWM)) (54) 
\Pz) = Y liO (0,4>) (55) 

where Yi rm (9,(f>) are spherical harmonics. From this we 
can construct the p xy orbitals as 

\p%) = -±= (\ Px ) ±i\ Py )) = TYi,±i(M) (56) 

which correspond to the two a + and er~ sub-bands, as we 
will call them in the following. If we introduce rotated 
spherical coordinates we can define the following set of 
transformations between the unrotated and the rotated 
frame 

sin 9' cos 4>' — sin 9 cos <fi (57) 
sin 9' sin <\> = cos 7 sin 9 sin tj> + sin 7 cos 9 (58) 
cos#' = — sin 7 sin 9 sin <j) + cos 7 cos 9 (59) 

where 7 is the misorientation angle (see Fig. 2J. Now 
we can calculate the overlap integrals for the different 
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orbitals as 

(Pz\Pz(l)) 
(pty\Pz{l)) 

(pty\Pxy(l)) 
(Pxy\pty(l)) 



cos 7 (60) 
T (Yl ±l {e,4>)Y 1 ,^'^'))n 

± 7I sin7 ^ 

(YZ ±1 (e,<t>)Y lt ± 1 {0',tf/)) a 

1(1 + cos 7) (62) 

-(Yi Tl {e^)Y^ ±l {e',<p')) n 

1(1 -cos 7 ) (63) 



where we have used the abbreviation — 

Jq d6 J Q 2lr d(j) sin 8 . . . . The missing integrals can be found 
by simple symmetry considerations. If we now calculate 
the total transmission probability for one particle for ex- 
ample from the er + -band into one of the other bands we 
get 

\{pty\pty{l))\ 2 + \{pty\Pxy(l))\ + \{pty\Pz{l))\ 2 

= -(1 + COS7) 2 + -(1 - cos 7 ) 2 + - sin 2 7 = 1 (64) 

as to be expected. We find the same total transmission 
propability for the other bands as well. If we multiply the 
overlap integrals with a global transmission factor t de- 
scribing the overall reflection and transmission properties 
of the point contact we can write the angular dependent 
effective hopping amplitudes as 



T rT + rT + 



-(1 + cos 7) 
77 (1 — cos 7) 



r n7T = t cos 7 
We can write them in matrix form as 



V2 



sin 7 



(65) 

(66) 

(67) 
(68) 

(69) 



The effective hopping amplitudes enter the transfer- 
matrix equation as described in section [n] In Fig. [21 
we show the angular dependence of the effective hopping 
amplitudes. Note, that interband transmission vanishes 
for 7 = 0, when there is no misalignment. On the other 
hand, for 7 = 7r/2 interband transmission becomes max- 
imal, however, intraband transmission from the tt band 
into the tt band r TT vanishes. This peculiar behavior is 
a direct consequence of the orthogonality of the a and tt 
bands. 
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FIG. 2: The effective hopping amplitudes between the differ- 
ent bands as functions of the misorientation angle a. Here we 

have added the contributions from the two different a bands 

1 1/2 



defining r CT(T = [2 1 
as before. 



= 2r CT 



and tv 



V. RESULTS 

In this section we will discuss the dependence of the 
local density of states, the critical current and the differ- 
ential conductance on the misorientation angle 7 between 
the two electrodes of a superconducting MgB 2 point con- 
tact. In the following we assume that the gap value in the 
vicinity of the point contact is not changed so it can be 
taken as the bulk value. All energies will be normalized 
to the larger gap value while we will take the ratio 
of the two gaps to be = ^A^^ First we will cal- 
culate the local quasiparticle density of states (LDOS) in 
the vicinity of the point contact. In equilibrium without 
an applied voltage we can calculate the LDOS in band a 
as 



N { F a) 
2iri 



I a,+ ^ ' ck,- 



(70) 

using the expressions of Eqs. 12">l) and (|26|) to calculate 
the quasiclassical Green's function on the incoming and 
the outgoing trajectory. In Fig.|21thc LDOS is shown for 
two different contacts. In the first case we assume that 
both electrodes are equally orientated (7 = 0) and we 
have no interband tunneling. Then the different bands 
remain decoupled and we find a standard s-wave spec- 
trum in the a bands as well as in the tt band. In the sec- 
ond case we assume a maximum misorientation (7 = §) 
and we find peaks at the position of the tt gap in the a 
band spectrum and vice versa due to interband hopping 
processes. In both cases the transparency parameter t 
was chosen to be t = 0.1. In Fig. 0] we present cal- 
culations of the temperature dependence of the critical 
current for the same parameters. For two equally ori- 
entated electrodes we find a very good agreement of the 
quasiclassical results for low transparencies with the an- 
alytical approximation given by Eq. Q35[l. This is clear 
since for 7 = we have a decoupled system and the crit- 
ical current can be understood as a weighted sum of two 
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FIG. 3: The quasiparticle density of states at a point contact 
calculated within a two band model. The full line shows the 
DOS of the quasiparticles in the a bands while the dashed 
line shows the DOS of the quasiparticles in the tt band. Both 
curves are normalized to the normal conducting DOS in the 
corresponding bands. 



single gap critical currents, that are well described by the 
Ambegaokar-Baratoff formula. For a maximum misori- 
entation of 7 = 5 we find a deviation from the classical 
behaviour and a nearly linear decrease of the critical cur- 
rent over a wide range of temperatures. In both cases we 
have taken the normalized density of states at the Fermi 
energy to be = jj and Np = i corresponding 

to band structure calculations for MgB2i22 The deviation 
from the Ambegaokar-Baratoff result seen in Fig. 0Jj can 
be attributed to interband transmission processes. In 
the next figures we consider a point contact with an ap- 
plied voltage. The current j and the differential conduc- 
tance g can be normalized using the conductance of the 
contact in the normal state gN ■ Employing Eq. i|28|) and 
Eq. (|36|l the normal state conductance for quasiparticles 
of the band a can be calculated as 



g N = 4n 2 eN { F a 



V 1 + 



2 2\ 
TT V 



1 



FIG. 4: Critical current as a function of temperature calcu- 
lated within the transfer-matrix approach for two different 
orientations of the electrodes (filled points). The dashed line 
gives the corresponding two-band Ambegaokar-Baratoff ap- 
proximation. 



the case of low transparency we can neglect multiple An- 
dreev reflections and use a simple model that can be de- 
duced from Eq. expanding it for small values of v 
and keeping only terms in second order. In this case only 
the f-matrix components f±i,o have to be calculated and 
they reduce to i±i,o = ^±i,o- Then all the different terms 
in Eq. (|36|1 can be summed up to 



(71)" 

and the total conductance is a sum over all band com- 
ponents gN — J2a9N- For high transparencies it shows 
a strong dependence on the misorientation angle 7 while 
it becomes independent of 7 in the limit of low trans- 
parency. In this limit also the calculation of the non- 
equilibrium properties can be considerably simplified. In 



J, 



■L = 4a 3 <<*.«>„ 



a, a' 

de 



(72) 



x / y Pi(e) P a R (e - eV) (f F (e) - f F (e - eV)) 



where we have used that in 



[<?f(e) 



equilibrium g^ 



K 



(e) 



tion f F (e) - uc^2T 
ten the quasiparticle density 



1 

" 2ixi 



R.a / \ 

9j ( e ) 



A,' 
9j 



9ji e )\ /i ? ( e ) w hh a fermionic distribution func- 
tanh 7^ ■ Furthermore we have writ- 
of states as pf(e) — 

' (e) . In Fig. [5] we have calculated 

the current voltage characteristics for three different mis- 
orientation angles 7 within this model. The dotted and 
dashed lines show the contributions of the different bands 
where we have added the formally equivalent contribu- 
tions from the two a bands. For two equally orientated 
crystals on both sides of the interface (7 = 0) we have 
no interband hopping amplitudes and therefore we find 
two distinct increases in the current stemming from the 
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FIG. 5: Current as function of applied voltage for an S-I- 
S junction calculated in a multiband model with parameters 
adapted to MgEb (e.g. A^ CT ' = 3A^) and low transparency. 
The angle 7 gives the relative orientation of the two crystals 
on both sides of the interface and determines the effective 
strength of interband and intraband tunneling. 
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FIG. 6: Differential conductance as function of applied volt- 
age for an S-I-S junction calculated in a multiband model 
with parameters adapted to MgB2 (e.g. A' CT ' = 3A^) and 
low transparency. The angle 7 gives the relative orientation of 
the two crystals on both sides of the interface and determines 
the effective strength of interband and intraband tunneling. 



a band contribution at V = 2A^ and from the 7r band 
contribution at V = 2AW. This can be seen as a sim- 
ple addition of two decoupled single band contributions 
with different gap sizes. For 7 = j we find three dis- 
tinct increases, the first one for V = 2A^\ the second 
one for V = AM + A^) and the last one for V = 2A< <t) 
since the intraband as well as the interband hopping am- 
plitudes become finite. For 7 = ^ we find a maximal 
increase in the tunneling current for V = A^ + A^ 
and a smaller contribution at V = 2A^ as expected 
from the angle dependence of the hopping amplitudes. 
Intraband processes between the it bands do not exist 



for this orientation since the ^ rotated p z orbitals are 
orthogonal. In Fig. [H] we show the corresponding differ- 
ential conductance as a function of voltage obtained by 
numerical differentiation of the current voltage charac- 
teristics in Fig. [3] Here the increase of the tunneling cur- 
rent due to tunneling between the high density of states 
at the gap edges results in distinct peaks in the differ- 
ential conductance. In comparison to Fig. [S] and Fig. 
the situation for point contacts with high transparency is 
totally different. Here the occurrence of Andreev reflec- 
tions dominates the current voltage characteristics. First 
we want to consider an N-I-S junction with multiband 
electrodes on both sides of the interface. This situation 
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FIG. 7: Differential conductance as function of applied volt- 
age for an N-I-S junction calculated in a multiband model 
with parameters adapted to MgB2 (e.g. A'' 7 ' = 3A' W ') and 
maximal transparency t = 1/tt. The angle 7 gives the relative 
orientation of the two crystals on both sides of the interface 
and determines the effective strength of interband and intra- 
band tunneling. 



FIG. 8: Differential conductance as function of applied volt- 
age for an S-I-S junction calculated in a multiband model 
with parameters adapted to MgB2 (e.g. A' CT ' = 3A^) and 
high transparency t — 0.1. The angle 7 gives the relative ori- 
entation of the two crystals on both sides of the interface and 
determines the effective strength of interband and intraband 
tunneling. 



could for example be experimentally realized by a re- 
duced critical temperature in one of the electrodes due 
to defects. Andreev reflection processes carry twice the 
charge as the corresponding one electron process. Thus 
we would expect in the case of maximum transparency 
an increase of the differential conductance for low volt- 
age leading to a value that is two times larger than the 
high voltage asymptotics. This case is realized in Fig. [7| 
for 7 = and the results are in good agreement with the 
BTK-results presented by Brinkman et alii&. If we as- 
sume a misorientation of the electrodes we find a decrease 
of the differential conductance for low energies that can 
be explained by a reduction of the current in the tt band 



as a result of an effectively smaller transparency of the 
barrier. 

In the last figure (Fig. [8} we show the differential con- 
ductance for an MgB2 S-I-S point contact calculated 
within the non-equilibrium quasiclassical theory. For 
7 = we find distinct conductance peaks at V — 2Aw 
and V = 2A^ and also the well known subgap struc- 
tures at V — 2A j ' i ) with n € N that appear due to mul- 
tiple Andreev reflections. If we rotate the electrodes we 
find a more complicated situation. For 7 ~ j the peak 
at V = 2A( <J ) shrinks visibly as expected due to a re- 
duced intraband hopping probability while the peak at 
V = |A( CT ) = 2AW increases. This can be easily un- 
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derstood since not only hopping processes between the 7r 
bands contribute to this peak but also multiple Andreev 
reflections corresponding to V = ^2A^ and interband 
processes with V = (A^ + AW)/2, For 7 = § the 
peak at V — 2A' <T ' has nearly vanished while we find 
at V = 2A( 7r ) a dramatic reduction of the differential 
conductance. Instead we find an increase of interband 
processes corresponding to V = A( CT ) + A^) = lA^J 
(arrows). This peak at AW +AW is a characteristic sig- 
nature of interband transmission. It is absent for 7 = 0, 
when the electrodes are aligned, and it becomes maxi- 
mum for 7 = 5' ^ n experimental observation of this 
effect would provide direct evidence for the orthogonal- 
ity of the wavefunctions belonging to the 7r and a bands. 

VI. CONCLUSIONS 

Due to its multiband character MgB2 Josephson junc- 
tion promise to show rich and interesting physical proper- 
ties. We presented calculations of current-voltage char- 
acteristics of both S-I-S and N-I-S contacts. We de- 
veloped a simple model to determine the effective intra- 
band and interband hopping amplitudes at the interface 
from a microscopic point of view. We used these model 
parameters within a multiband transfer-matrix descrip- 
tion based on the work of Cuevas and Fogelstrom and 
Kopu et al. to study the dependence of the tunneling 
currents on the relative orientation of the two electrodes. 
Especially due to the effect of multiple intraband and in- 
terband Andreev reflections we were able to identify dra- 
matic changes in the differential conductance for different 
electrode orientation. For N-I-S junctions with equally 
orientated electrodes our result agrees with the BTK- 
results of Brinkman et al. For S-I-S junctions we found 
subgap structures at different voltages corresponding to 
different combinations of the two gaps. The features sen- 
sitively depend on the misorientation angle of the junc- 
tion. We hope that especially the angular dependence of 
the differential conductance can be experimentally ver- 
ified for point contacts on MgB2 single crystals which 
would provide direct evidence for the orthogonality of 
the two bands. 



formed into an expression that only contains the retarded 
and advanced components of the off-diagonal transfer 
matrix elements Ilr and i R L. 



1. Relation between the Keldysh and the retarded 
and advanced components 

First we show that the Keldysh component of the T- 
matrix can be replaced by the corresponding advanced 
and retarded components. In the following we will show 
the relations for the Gor'kov Green's functions but of 
course they are also valid for the components of the quasi- 
classical Green's propagator. Starting with Eq. Q, writ- 
ing it out in Keldysh-space we get for the retarded and 
the Keldysh components of the T-matrix the following 
relations 

f R = d + i,oC*of R (Ai) 

and 

t K = ~vo^ot K + ^oi A ) (A2) 

Solving the first equation for 1~ R and the second equation 
for T K we get the following expressions 

f R = (l-io^Y'oi (A3) 

and also 

t K ^(i-ioC^y 1 0*0^0^ (A4) 

Substituting the first equation into the second equation 
we end up with the important relation 

f K = f R o o f A (A5) 

With this equation we are able to replace the Keldysh 
component of the transfer matrix by the Keldysh com- 
ponent of the decoupled Green's function, that can be 
easily calculated within the quasiclassical theory. 
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APPENDIX A: CALCULATION OF THE 
CURRENT EXPRESSION 

In this appendix we will show how the commutator 
formula for the spectral current density can be trans- 



2. Relations between the reflection and the 
transmission components of the transfer-matrix 

In a second step we will replace the reflection compo- 
nent T^L by the corresponding transmission components 
Tlr and Trl. Comparing the expressions for the trans- 
mission and reflection components we find 

TLL = (I - V LR o Qoo^r o Vr L ° Goo, l) 

°VLR ° QocM ° v RL (A6) 
Tlr = (l - v LR o Qoom o v RL o Qoo^l) 1 O V LR 
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and therefore 

Tll = T LR o o v RL (A7) 

To calculate the commutator appearing in the current 
formula we are interested in expressions containing Tll o 
Goo.l and Qoo.l o I'll. To find them we start with 

= v LR o Qoo^r o v RL o(l + l o T LL ) (A8) 

and performing 

Tlx o(l + 0^ o T LL ) 1 = v LR o g^R o v RL (A9) 

we are now able to make use of the very important rela- 
tion ao (l + Joa) = (1 + a o b) o a and we get 

(l + I'll 3oo,l) 1 = i> LR og oo . R ov RL (A10) 

or 

TLL = (l + "TLL O Ooo.l) ° t'Lfl o ^00, _R ° "RZ, (All) 

Now we are able to replace the left hand side of Eq. (|A7|) 
with the expression above and after removing g ociy rovrl 
we find 

f LR = (1 + f LL o C?^) o t, ii? (A12) 

With basically the same idea we can find a relation be- 
tween Tll and Trl- We start with 

TRL = (l - V RL o g^x o v LR o Soo.i?) °v RL (A13) 

Multiplication of ^00. l from the right and using again 
the relation 00 (l + 6oa) 1 = (1 + a o 6) 1 o a leads to 

Tki Coo,L — VRL ° <?oo,L 

(I - v LR o 0^ o v RL o Coo^l) 1 (A14) 

Another multplication from the right, now with Vlr, let 
us compare the right hand side of the equation with the 
definition of 1~lr and we find 

TrL ° £oc,L ° V L R = 1> RL ° £00, L ° T LR (A15) 
TLR OocM ° V RL = V LR O 0oo,R o Trl (A16) 

where the second equation has been found by interchang- 
ing L and R. With this connection between Tlr and Trx 
we can finally find another useful relation by substituting 
the expression for Tlr from Eq. (|A12|I on the right hand 
side and removing the 600, l v lr factor 

tiu.=*RL° + °fz.i) (A17) 
3. Calculation of the spectral current density 

The spectral current density can be calculated by 
Eq. H29|) . In the following we will use the abbreviation: 

K* = [f L L,Gao,L]* (A18) 



If we evaluate the commutator and perform the matrix 
multiplications in Keldysh space we get 

K? = T'll o QLl + Tll o Ci, L (A19) 

In the next step we can replace the Keldysh components 
of the transfer-matrices by its corresponding retarded 
and advanced components employing the relation of ap- 
pendix IA II 

*LL = 7L" ° G*,L ° t£ + ^LR O <?£,H O 7& (A20) 

Now the commutator takes the form 

*f = t L \og^ L o{l + t£og^ L ) 
-(i + g^ L ot L l)og^ L ot L A L 

+Tl R r°G£,r°Trl°G£,l (A21) 

In the last step we will replace the reflection components 
of the transfer matrix by the corresponding transmission 
components using the relations of appendix IA 21 We start 
by substituting the brackets using Eqs. I|A12(I and (|A17(I 



7m *Zi = (l + 7& ° #U) (A22) 

and 

««i ° - (1 + g«, L ° t£) (A23) 

The remaining two reflection components can be elimi- 
nated employing Eq. (|A7|) and Eq. i|A16p leading to 

= T L A R ogA R ov RL , f L \ = v LR og* R of£ L (A24) 

and we can finally write down the commutator corre- 
sponding to the spectral current density as 

K? = VLRog^ R of* L ogK L of L A R ov- R 

~ v rl Trl ° g*, L ° Tlr o g A jR o v RL 
4VC ot A°^ (A25) 

Performing the quasiclassical limit this expression defines 
the spectral current as it is given by Eq. 
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